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ABSTRACT 

Aims. To incorporate background subtraction into the Bayesian Blocks algorithm so that transient events can be timed accurately and 
precisely even in the presence of a substantial, rapidly variable, background. 

Methods. We developed several modifications to the algorithm and tested them on a simulated XMM-Newton observation of a bursting 
and eclipsing object. 

Results. We found that bursts can be found to good precision for almost all background subtraction methods, but eclipse ingresses 
and egresses present problems for most methods. We found one method that recovered these events with precision comparable to 
the interval between individual photons, in which both source and background region photons are combined into a single list and 
weighted according to the exposure area. We have also found that adjusting the Bayesian Blocks change points nearer to blocks with 
higher count rate removes a systematic bias towards blocks of low count rate. 

Key words, methods: numerical - (stars: ) binaries: eclipsing - X-rays: binaries 


1. INTRODUCTION 

The detection and precise timing of transient events is of im¬ 
portance in all areas of astrophysics. In some applications, for 
instance X-ray astronomy, the desired timing accuracy is com¬ 
parable to the interval between the arrival times of individual 
photons. In this situation it is necessary to identify unambigu¬ 
ously which photons mark the end of the pre-transient phase and 
the beginning of the post-transient phase and, if possible, esti¬ 
mate without systematic bias where between those photons the 
change occurs. 

One method t hat has attracted recent interest is t he Bayesian 
Blocks algorithm (IScarglelll998l : IScargle et al.ll2013ll . This algo¬ 
rithm may take as its input a list of photon arrival times, such as 
those produced by X-ray telescopes. The algorithm can also han¬ 
dle other data modes, but in this work we consider exclusively 
time tagged event data. The observation is then divided into cells, 
intervals containing a single photon, and periods of time of nom¬ 
inally constant count rate are produced by finding an optimum 
number and placement of blocks of consecutive cells. Given a 
prior on the number of blocks, the algorithm finds an objectively 
optimal set of blocks describing the observation, and within each 
block the arrival times of the photons are consistent with a Pois¬ 
son process with constant rate. The Bayesian Blocks algorithm 
is ideal for finding and accurately timing transients, as a change 
point between blocks, typically placed half-way between the last 
event of one block and the first event of the next, will be found 
if and only if the rate of arrival of photons detectably changes. 
Transient timing can therefore, in principle, be performed with 
a precision comparable to the interval between individual pho¬ 
tons, and clearly no better precision is possible. This property 
makes the Bayesian Blocks algorithm an attractive candidate for 
transient timing. 

Other binning methods exist which, similarly to the Bayesian 
Blocks method, place bin edges at locations o nly where the data 
justifies it (e.g.. lKnuthll2006c lBelangerll2013l) . but we have not 


considered them in this paper. See also iBurgessI (l2014t) for a 
comparison of binning methods, in the context of gamma ray 
burst timing. 

The motivation for this paper is accurate and precise eclipse 
timings of cataclysmic variables (CVs) and low-mass X-ray bi¬ 
naries (LMXB) using data from the XMM-Newton X-ray obser¬ 
vatory. Much of the paper is presented with these applications in 
mind, but the methods and conclusions developed are not spe¬ 
cific to this astronomical field and are not even limited to astron¬ 
omy. We consider two effects that may affect accurate timing 
measurements using the Bayesian Blocks algorithm. 

Firstly, is the location of the cell edges between photons 
optimal? While most applic ations put the cell edges half way 
between tw o pho tons (e.g., IScar ^Il998t IScargle et al.ll2013l 
llvezic et^l2014h . there is no reason why the edge cannot be 
placed at any point between the photons. In ^we describe a sys¬ 
tematic bias affecting the “half way” placement and suggest an 
adjustment that removes this bias almost entirely. We also con¬ 
sider the case of a continuously varying source intensity, such as 
those caused by an extended emitting spot moving into or out of 
eclipse. 

Secondly, XMM-Newton observations are frequently affected 
by soft proton flares originating in Earth’s magnetosphere that 
contribute a substa ntial, and often rapidly variable, background 
(iLumb et al.ll2002h that needs to be removed to recover the true 
variability of the source. Our data therefore consists of photon 
lists: one extracted from a region of sky surrounding the source 
and containing both source and background photons, and one 
taken from a source-free region of sky containing background 
only. It is not obvious how to perform background subtrac¬ 
tion using the Bayesian Blocks algorithm. Removing individual 
source region p hotons according to the background rate has been 
suggested (e.g.. lStelzer et al.|l2007h . but we would prefer not to 
discard any data. An attempt at weighting individual photons in 
a manner very similar to the ones de veloped in this paper has 
been recently applied successfully bv iMossoux et al.l (120151) . In 
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^we generate theoretical source and background light curves, 
containing transient events occurring at known times. In ^we 
propose several modifications to the Bayesian Blocks algorithm 
to allow for background subtraction, and evaluate them for their 
ability to recover the transients in the model light curve in the 
presence of a strong, rapidly varying, background. 


2. ANALYSIS OF ERRORS 


In this paper we use the geometric prior suggested by 
IScargle et al.l (l2013h . This prior contains one adjustable param¬ 
eter po, reflecting the balance between suppressing spurious 
change points and retaining genuine ones. Throughout this paper 
we have taken po - 0.01, which selects against false positives at 
the 99% level. 

Most implementations of the Bayesian Blocks algorithm 
place the cell bo undaries at the midpoint between two succes¬ 
sive e vents (e.g.. lScarglelll998HScargle et al.ll2013Ulvezic et al.l 
|2014|) . This will bias the location of change points slightly to¬ 
wards blocks of low count rate. Consider a Poisson process that 
has rate ro for times f < 0 and rate ri for f > 0, with ri < tq. 
Suppose the last event of the first segment and the first event of 
the second segment occur at times to and fj respectively. Since 
they are both Poisson processes, to and fi will have mean val¬ 
ues of -l/2ro and l/2ri, and the Bayesian Blocks change point, 
since it falls halfway between the two, will have mean value 
(fQ - ^i)/(2r‘ori). This is positive, so the location of the change 
point is biased toward the block with the lower count rate. Figure 
[T] illustrates this bias. For applications such as eclipse timings 
in the study of cataclysmic variables, this bias can be a severe 
drawback, as the count rate during the eclipse is likely to be very 
low. The effect of this bias will be to cause the eclipse to appear 
shorter than it really is. Furthermore, if the pre- and post-eclipse 
blocks are different in count rate, the eclipse midpoint may be 
shifted. 

It is possible to improve the positioning of the change points, 
as follows: suppose the Bayesian Blocks algorithm has found a 
block of high count rate beginning at tg, containing no events, 
concluding with an event recorded at time to, followed by a block 
of lower count rate. We want to find the optimum location of the 
change point tcp. Since we are assuming a Poisson process, to is 
exponentially distributed with standard deviation 1 /ro where ro 
is the count rate in the block. It follows that, on average. 


1 

^cp — ^0 "Z ■ 

2ro 

We also have 
no 

ro - - 

^cp 


( 1 ) 

( 2 ) 


and these can be solved simultaneously to give 


^ _ 2«ofo - ts 

2no-l • 


( 3 ) 


The case of a rising count rate is very similar. 

We tested this new adjustment with a simulated event se¬ 
ries consisting of 100 events with count rate 3, up to f = 0, fol¬ 
lowed by another 100 events with count rate 2/3, and found the 
change point according to both the “half way” method and the 
one given by Equation |3] We repeated this test for 50,000 re¬ 
alizations of the series. The count rates here are dimensionless 
because the unmodified Bayesian Blocks algorithm is invariant 



Fig. 1. Illustration of the change point bias toward blocks with lower 
count rate. The two count rates, 3 and 2/3, are shown as the blue line, 
with events placed according to a Poisson process shown as circles. 
The half-way point between the two terminal photons is shown as a 
solid vertical line, as the change point placed according to Equation |3 
is shown as a vertical dotted line. It is clear that the “half-way” change 
point method is biased towards blocks with lower count rate but that 
the adjusted change point location gives a better estimate of the true 
location of the change point. 


under changes in the unit of time. The mean location of the re¬ 
covered change point using the “halfway” change point location 
was tcp - 1.36 + 1.83, comparable to the average spacing of pho¬ 
tons in the lower count rate segment and indicating a substantial 
bias toward the lower count rate. For the new method, the mean 
location was tcp - 0.461 + 1.58, an improvement by a factor 
of three. The small bias that remains is entirely due to the algo¬ 
rithm occasionally mistaking the first event of the second seg¬ 
ment for an event in the first segment, if it happens to be placed 
close to the true change in count rate. The converse, mistaking 
an event in the higher count rate segment for one in the lower 
count rate segment, is also possible but much less likely since 
this requires several photons in the higher count rate segment to 
be badly placed. 

We repeated the experiment, discarding all trials with 
misidentified photons, and found that the mean change point lo¬ 
cations for the “half way” and adjusted change point locations 
were 0.66 + 0.45 and -0.044 + 0.20 respectively. The adjusted 
change point method suffers from essentially no bias, except for 
the uncertainty in determining which photons belong to which 
blocks. 

Adjusting the location of change points by this method will 
only be useful if the uncertainty in their location is small com¬ 
pared to the reduction in bias. To investigate under what cir¬ 
cumstances the adjustment is useful we performed the follow¬ 
ing tests. We produced many segmented light curves as above, 
consisting of 100 events with count rate ri = 1, followed by 
100 events of count rate ro between 0.001 and 1. The division 
between the two count rates is at f = 0. We produced 1,000 
realizations of each of these light curves and performed the 
Bayesian Blocks algorithm on them, with and without the new 
change point adjustment, repeating trials where no change point 
was found at all. Then we compared the improvement in bias 
achieved (i.e., the difference of their absolute magnitudes) and 
compared it to the remaining uncertainty. The results are shown 
in Figure 12] 
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Fig. 2. Reduction in bias by employing the change point adjustment 
(black points), and remaining uncertainty in change point location (grey 
points), plotted against the lower of the two count rates. The improve¬ 
ment is larger than the remaining uncertainty for to <0.15 


The change point adjustment provides an improvement in 
bias greater than the remaining uncertainty for ro smaller than 
about 0.15 but no useful improvement for higher tq. This be¬ 
haviour can be understood by considering the sources of uncer¬ 
tainty and bias for the two methods. For the unadjusted method, 
both bias and uncertainty are dominated by the placement of in¬ 
dividual photons in the lower count rate segment. A misidenti- 
fied photon contributes relatively little to the bias because, to be 
misidentified, it must be placed within about 1 /ri of the change 
in count rate, and this is smaller than I/tq. For the adjusted 
method the bias consists almost entirely of misidentified pho¬ 
tons, and the uncertainty is partially due to misidentified photons 
and partially due to the placement of individual photons in the 
ri segment. It follows that the bias and uncertainty in the un¬ 
adjusted method decreases with increasing ro, because 1/ro de¬ 
creases, and both the bias and uncertainty in the adjusted method 
increase because misidentified photons become more and more 
likely. 

The adjustment of change points according to Equation[3]is 
therefore most useful when one of the blocks is less than 15% as 
intense as the other. For the remainder of the paper, we will be 
using this method to determine the location of the change points. 

2.1. Continuously varying count rate 

In studies of eclipsing CVs, it is useful to determine or constrain 
the size of the emitting accretion spot on the white dwarf. As the 
spot passes behind the secondary star, its observed flux declines 
steadily to zero. If this decline is gradual enough, the Bayesian 
Blocks algorithm will produce blocks of intermediate count rate 
during the ingress or egress. If the decline is too rapid, the al¬ 
gorithm will not resolve it and produce only an instantaneous 
change in count rate; the failure can still be used to constrain the 
size of the emitting region. 

We investigated this issue by simulating many hypothetical 
eclipse egresses, with a very low eclipse count rate of 10“^ and 
post-eclipse count rates R between 1 and 15. The egress itself 
was modelled as a linearly rising count rate with duration At be¬ 
tween 0 and 50 and the fiducial time t - 0 was placed at the 
midpoint of the egress. For each set of R and At we made 1,000 
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Fig. 3. Probability of resolving an eclipse egress with post-eclipse count 
rate R and eclipse egress duration At. Contours are shown at the 50%, 
90%, and 99% significance level (solid, dotted, and dashed lines respec¬ 
tively). The probability of finding more than one change point during 
egress increases with increasing duration and post-eclipse count rate. 

realizations and found the Bayesian Blocks change points. Fig¬ 
ure [3 shows the probability that the egress will be resolved, that 

is, that the Bayes Blocks algorithm finds more than one change 
point in the egress. It is clear that this probability increases with 
increasing egress length, and with increasing post-eclipse count 
rate, as one would expect. 

The usefulness of this analysis to eclipse timing is clear. If 
the ingress or egress is not resolved, one can place constraints on 
the size of the emitting spot by finding the smallest At for which 
that ingress or egress would have been resolved. For instance, 
if an eclipse egress in a real observation is not resolved and the 
post-eclipse count rate is 6 photons/s, it can be seen from Figure 
[3 that the egress has a shorter duration than 30 s with approxi¬ 
mately 99% confidence. 

For simulations where the egress was not resolved, that is, 
where only one change point was found for the egress, we found 
the mean position of that change point. The results are shown in 
Figure m Three behaviours are evident in this Figure. For short 
eclipse egresses, the situation closely resembles an instantaneous 
jump in count rate and the bias on the position of the change 
point is accordingly scattered around zero, with the scatter ap¬ 
pearing large on the Figure because At is small. For intermedi¬ 
ate durations, where the duration is long compared to the photon 
separation in the post-eclipse block, the change point is triggered 
earlier than the half-way point of the egress. There is a bias of 
about 25 to 33%. Finally, when the post-eclipse count rate is high 
and the egress is long, then the algorithm never fails to resolve 

it, as indicated by the absence of such points in Figure|4] 


3. SIMULATED DATA FOR BACKGROUND 
SUBTRACTION 

In this section we devise and apply several methods for perform¬ 
ing background subtraction using a Bayesian Blocks approach. 
We have tested the various methods in this paper on a hypo¬ 
thetical XMM-Newton observation of a binary system exhibiting 
transients in the form of eclipses and bursts. The observation of 
the source is affected by a soft proton flare, contributing a sub¬ 
stantial background of photons that increases in magnitude over 
the observation, and is highly time-variable. 
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Fig. 4. Bias in the eclipse egress change point location as a fraction of 
egress duration, for simulations in which the egress was not resolved. 
Colors indicate the post-eclipse count rate. 


We generated two lists of photon arrival times. The first is 
assumed to be taken from a region surrounding the source, con¬ 
taining photons from the source itself and from the soft proton 
flare superimposed upon it. The second list contains only flare 
photons, taken from a region 4.123 times larger than the source 
extraction region. This number was arbitrarily chosen. Thus, 
photons from the background extraction region will be given a 
weight of Wb - 1/4.123 compared to source region photons. 
We will attempt to recover the source light curve by subtract¬ 
ing the soft proton flare from the total light curve. The ability 
of our background subtraction methods to recover the timings of 
the transient events is our measure of the suitability of the back¬ 
ground subtraction methods. We summarize the results in Table 
|T]and Figure [13] 

The source and flare light curves are described in ^3.11 and 
^3.2l resr)ectivelv. 


3.1. Source light curve 

We have taken the behaviour of the LMXB EXO 0748-676 as 
a model for our simulated source light curve. This system ex¬ 
hibits eclipses of 8.3 minute dur ation, with an orbital period 
of 3.82 hours (iParmar et al.lll98^ . during which the X-ray flux 
from the source drops to zero. It also undergoes type-I X-ray 
bursts, during which the X-ray flux incr eases almost instanta - 
neously to many times its pre-burst level dGottwald et al.l[T98^ . 
When neither eclipses nor bursts are present, the source emits 
with a count rate of abou t 3 counts/s in XMM-Newton observa¬ 
tions dHoman et al.ll200^ . 



Time (s) 


Fig. 5. Total source light curve (top panel), background (middle panel), 
and zero-background source light curve (bottom panel). The high inten¬ 
sity of the background signal is due to the higher exposure area (4.123 
times greater than the source area). 


3.2. Background light curve 

The background consists of three segments. First is a linear rise 
from a count rate of 0 to 24, to f = 7500, followed by a quadratic 
fall back to zero. These two segments are intended to investigate 
whether the presence of a rising or falling background contribu¬ 
tion biases the timings of transients. Beginning at 15,000 s is an 
oscillatory signal with increasing amplitude and frequency, de¬ 
signed to approximate the behaviour of the soft proton flares that 
affect observa tions by X-ray inst ruments such as XMM-Newton 
and Chandra dLumb et al.ll2002h . Its count rate is given by the 
formula 


R(t) = 2.0 X lO^^f^ 


1 + sin 


1 


5000 


1.25 


log, h 


photons/s, (4) 


The simulated source light curve consists of a constant 
persistent intensity of 3 counts/s. Superimposed upon this are 
eclipses of duration 498 s, recurring with a faster orbital period 
of 3,035 s, and bursts recuri'ing with a period of 2545.6 s. The 
bursts have peak intensity of 30 counts/s and exponential decay 
time of 24 s. The transient events in the light curve are therefore 
(see Figure |5]l ten bursts, ten eclipse ingresses, and ten eclipse 
egresses. All transients occur instantaneously, and two bursts are 
“missing” due to falling inside an eclipse. The times of these 
events are given in Table [1] There are 154,992 photons in this 
series. 


where ta-t- 15000. There are 300,725 photons in this series. 

4. BACKGROUND SUBTRACTION 

4.1. Constant Cadence 

Before analyzing various implementations of the Bayesian 
Blocks algorithm, we study a case in which the locations of 
the change points are not optimally determined by the data, as a 
comparison. The simplest method is to simply group the source 
and background photons into equally spaced bins whose width 
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Fig. 6. Source, background, and background subtracted light curves 
(top, middle and bottom panels respectively) accumulated into bins of 
predetermined width (10 s) and position. This method is effective at re¬ 
moving the background but its ability to time transient events is limited 
by the previously set bin size. 
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Fig. 7. Direct subtraction of the Bayesian Blocks representation of the 
background from that of the source. This method produces one change 
point for each change point in the two Bayesian Blocks representations 
from which it is derived. 


and location are fixed beforehand, then subtracting the latter 
from the former to obtain the background subtracted light curve. 
In Figure |6] we show the two photon series, and their difference, 
counted in bins of 10 s width. 

A property of the constant cadence method is that signals 
with more rapid variability than the cadence will tend to be aver¬ 
aged out. This property makes it useful for subtracting a rapidly 
varying background. However, fast transients like burst rises and 
eclipse in- and egress (instantaneous in our simulated data) will 
not be timed accurately unless they fortuitously happen to fall on 
or near a bin boundary. An instantaneous transient event falling 
well inside a bin will not appear instantaneous because it will 
produce one bin of intermediate count rate, giving the false im¬ 
pression of a gradual change; this phenomenon is clearly seen in 
panel b of Figure [121 In this regime, the uncertainty in the tim¬ 
ing of the transient event will be less than Af/2, where At is the 
width of the bins. 

A further disadvantage of the constant cadence method is 
that it is necessary to select the width of the bins beforehand. 
Vari ous estimates for the optim um bin width have been proposed 
(see iBirge & Rozenhoi3l2006l for a discussion), but they often 
suggest too f ew bins t o perfo rm timing analyses. The commonly 
used rule of ISturge^ (Il926h . taking 1 + log 2 V intervals, gives 
only 18 or 19 bins fo r the 154,992 source photons. Sim ilarly, the 
rules of lScottl (Il979t) and iFreedman & DiaconisI (Il981h give bin 
widths of hundreds of seconds. 


4.2. Bayesian Blocks- Direct Subtraction 

The most obvious way to perform background subtraction is to 
generate Bayesian Blocks light curves for the source and back¬ 
ground regions and subtracting the latter from the former, after 
normalizing the background series to account for the larger size 
of the extraction region. The results of this test are shown in 
Figure [T] and panel c of Figure [12] The resulting background- 
corrected light curve will contain a block boundary for every 
change point in both of the source and background Bayesian 
Blocks representations, and many of these will be spurious. 
Subtracting the background Bayesian Blocks light curve from 
the source light curve succeeds in removing much of the back¬ 
ground. However, the resulting light curve is very noisy, partic¬ 
ularly towards the end of the observation. 

IScargle et al.l (1201 3h provide a method for combining mul¬ 
tiple data series in such a way that the change points in each 
series line up. Such a process is useful for concurrent observa¬ 
tions of the same transients by different instruments, and will 
not produce superfluous change points when they are added or 
subtracted, but is not suitable for data series such as ours, which 
contain completely different transients. 

If we do not wish to produce a Bayesian Blocks repre¬ 
sentation of the background, the background can be subtracted 
from the source Bayesian Blocks representation on the level of 
blocks, cells, or individual photons. These three approaches are 
described in ^4.3114.41 and l4.5l 
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Fig. 8. Background corrected light curve using the “weighted blocks” 
method. The background subtraction appears to be very effective, pro¬ 
ducing almost no noise even at the end of the observation. The eighth 
and ninth eclipses appear to have shallowly sloping edges, probably due 
to the presence of pulses in the background. 


Fig. 9. Background corrected light curve using the “weighted cells” 
method. There is considerable noise near the end of the observation, 
but all the transients appear sharply defined. 

cell contains not one photon, but 1 - hbWb photons, where is 
the number of background photons occurring within that cell. 
As the fitness functio n of each potential block (Equation 19 
of IScargle et al.ll2013h involves the logarithm of the counts in 
blocks, we must find a way of dealing with cells and blocks with 
nonpositive count rates. If the count rate is not positive we set 
the effective count rate in the logarithm to a small positive num¬ 
ber Smin (we have taken Smin = 1.0 x 10“^). The cells themselves 
are defined by the source photons, and therefore each contains 
at least one positively weighted photon, hopefully reducing the 
number of times we need to resort to this arbitrary countermea¬ 
sure. If a block edge separates two blocks of apparently negative 
count rate, we do not adjust its location according to Equation |3] 
because it is unclear what the spacing between two hypothetical 
photons in such a block should be. 

The results of this trial are shown in Eigurej^and panel e of 
Eigure[l2] The “weighted cells” method performs well through¬ 
out most of the observation, but towards the end it produces 
many short spurious blocks with implausibly high count rates. 
It also badly misses some eclipse timings (see panel e of Eigure 
[T2l and Table [TJ. 


4.3. Bayesian Blocks- Weighted Blocks 

This method is suggested by the observation that photons from 
the background region are more numerous, but individually carry 
little weight, compared to source photons. Here we find the 
Bayesian Blocks change points of the source photon list, and 
subtract from the source photon count rate the exposure area 
weighted count rate of the background photons falling into that 
block. That is. 


where CR is the background subtracted count rate, ns and hb are 
the numbers of source region and background region photons 
respectively, and L is the length of the block. The results are 
shown in Eigure [8] and in panel d of Eigure [12] 

Although this method does not produce as many spurious 
change points as the previous one, it is clear that it cannot give 
the locations of the transients as accurately, because it contains 
only the change points of the source series and these are poten¬ 
tially offset from their true locations by the variable background, 
as can clearly be seen in Eigure [T2l 

4.4. Bayesian Blocks- Weighted Cells 

In this approach we subtract the weighted background photons 
from each cell before the change points are found. Thus, each 


4.5. Bayesian Blocks- Weighted Photons 

This approach is similar to the previous one, except that we no 
longer define the cell edges by source photons only. Instead, we 
combine both photon lists into one, and each cell has a weight 
of either 1 or -Wb, depending on whether it contains a source 
or a background photon. As in the previous method, we do not 
adjust the location of the change point between two blocks with 
negative count rates. 

The problem with negative count rates is more prominent 
now, because there are many cells that contain a background 
photon and therefore have negative count rates. However, the 
placement of change points is potentially finer. The results are 
shown in Eigure [T0| and panel / of Eigure [12] There are many 
short blocks with absurdly high or low count rates near the end 
of the simulated observation, and even a few nonsense blocks 
near the beginning. The timing properties of this approach are 
excellent nonetheless. 

4.6. Bayesian Blocks- Iterated Bayesian Blocks 

This approach is the one dev eloped and successful ly applied to 
X-ray activity of Sgr A* by iMossoux et al.l (l2015h . It is effec¬ 
tively a hybrid of the weighted cells and direct subtraction ap¬ 
proaches. Here, a Bayesian Blocks representation is produced 
for the background and source light curves to obtain count rates 
for the source and background regions. The count rates are then 


Article number, page 6 of[Ii 



























































































H. Worpel and A. D. Schwope: Background subtraction and transient timing with Bayesian Blocks 



0 5000 10000 15000 20000 25000 30000 


Fig. 10. Background corrected light curve using the “weighted photons” 
method. Noise, in the form of short blocks with implausibly high or low 
count rates, is more prominent in this method, but the transients are 
sharply defined. 



Fig. 11. The method described in iMossoux et alJ ll201.'5h . The top panel 
shows the original photon weighting and the bottom panel shows the 
alternative photon weighting. The original weighting does not actu¬ 
ally subtract the background, but only attempts to locate the transients, 
which explains why the eclipses do not have count rates near zero. Both 
methods are only noisy near the end of the observation. 


used to calculate appropriate weights for photons in the source 
region and the algorithm is then run a third time. A photon falling 
within a block of source region count rate Cs and background 
region Cb is given a weight of m = CsliCs + Cb), where Cb 
is positive and corrected for exposure area. It is clear to see that 
weighting the photons in this way does not actually subtract the 
background, but is intended only to find the locations of change 
points. For this reason we also test this recipe with an alternate 
weighting, w - \ - CbICs- The results are shown in Figure [TT] 
and panels g and h of Figure [121 


4.7. Invariance to time units 

The unmodified Bayesian Blocks algorithm produces the same 
change points regardless of the choice of time unit. We have in¬ 
vestigated whether this desirable behaviour is preserved for the 
various modifications described above. We took the same pho¬ 
ton series and divided the arrival times by 3600 to express the 
arrival times in hours rather than seconds. Naturally the constant 
cadence, direct subtraction, and weighted blocks methods are 



23000 23500 23000 23500 


Time (s) 

Fig. 12. The region around eclipse #8 and burst #9. a) true source light 
curve (red) and background (blue), b) constant cadence, c) direct sub¬ 
traction, d) weighted blocks, e) weighted cells, f) weighted photons, g) 
iterated Bayesian Blocks with original weighting, h) iterated Bayesian 
Blocks with alternative weighting. The true locations of the transients 
are indicated by dotted vertical lines. All other methods locate the burst 
accurately but many fail to accurately time the eclipse ingress, which 
coincides with a pulse of the variable background. The “weighted pho¬ 
tons” approach performs very well in locating all transients. 


unchanged under this transformation, but all other methods pro¬ 
duced many more change points than previously. The “weighted 
photons” method, for example, produced 16,080 change points 
for the time-scaled photon series, compared to 971 for the orig¬ 
inal series. This behaviour is probably due to our safeguard 
against negative block count rates. Since we have achieved toler¬ 
ably good results with count rates of 1 to 10, but higher effective 
count rates produce very many spurious change points, it may 
be generally desirable to scale the time unit. We have also found 
that the number of change points produced is fairly insensitive 
to the choice of imin- 
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Transient 

Time (s) 

<1 

<1 

<1 

Ate 

Atf 

Atg 

Ath 

I-l 

1751.0 

-1.0 

0.4 

0.4 

0.1 

0.5 

1.7 

1.7 

1-2 

4786.9 

3.1 

-0.2 

-0.2 

0.4 

0.1 

0.4 

0.4 

1-3 

7822.7 

-2.7 

4.0 

4.0 

-3.0 

3.7 

4.4 

4.4 

1-4 

10858.6 

1.4 

-0.1 

-0.1 

-2.3 

-0.1 

0.7 

0.7 

1-5 

13894.5 

-4.5 

0.1 

0.1 

-0.7 

0.1 

8.9 

8.9 

1-6 

16930.3 

-0.3 

0.1 

0.1 

-3.2 

0.1 

2.4 

2.4 

1-7 

19966.2 

3.8 

-0.2 

-0.2 

-0.2 

0.0 

6.7 

6.7 

1-8 

23002.1 

-2.1 

-8.3 

68.3 

-3.4 

0.5 

68.9 

-8.2 

1-9 

26037.9 

2.1 

-5.0 

25.2 

0.8 

0.5 

25.0 

5.5 

I-IO 

29073.8 

-3.8 

1.1 

1.1 

-0.6 

1.2 

2.1 

2.1 

E-1 

2249.0 

1.0 

-0.0 

-0.0 

0.5 

0.2 

-1.4 

-1.4 

E-2 

5284.9 

-4.9 

0.3 

0.3 

0.5 

-0.3 

-0.2 

-0.2 

E-3 

8320.7 

-0.7 

-1.3 

-1.3 

2.0 

0.0 

-1.8 

-1.8 

E-4 

11356.6 

3.4 

-0.3 

-0.3 

0.1 

-0.0 

-3.5 

-3.5 

E-5 

14392.5 

-2.5 

-0.1 

-0.1 

0.7 

0.9 

-23.2 

-23.2 

E-6 

17428.3 

1.7 

-1.4 

-1.4 

0.4 

0.1 

-2.5 

-0.5 

E-7 

20464.2 

-4.2 

-4.9 

-4.9 

-0.1 

-4.4 

-5.4 

-5.4 

E-8 

23500.1 

-0.1 

-16.4 

-16.4 

40.9 

-1.1 

-16.8 

-14.1 

E-9 

26535.9 

4.1 

16.2 

44.0 

35.2 

-0.3 

44.4 

16.0 

E-10 

29571.8 

-1.8 

-0.3 

10.7 

15.5 

-0.7 

21.5 

-0.3 

B-1 

800.0 

0.0 

-0.0 

-0.0 

0.0 

0.1 

-0.3 

-0.3 

B-2 

3345.6 

4.4 

-0.0 

-0.0 

0.1 

0.1 

-0.2 

-0.2 

B-3 

5891.2 

-1.2 

-0.0 

-0.0 

0.1 

0.0 

-0.3 

-0.3 

B-4 

8436.8 

3.2 

-0.1 

-0.1 

-0.6 

-0.6 

-0.2 

-0.2 

B-5 

13527.9 

2.1 

-0.1 

-0.1 

0.0 

0.0 

-0.4 

-0.4 

B-6 

16073.5 

-3.5 

0.0 

0.0 

0.2 

0.2 

-0.2 

-0.2 

B-7 

18619.1 

0.9 

-0.0 

-0.0 

0.0 

0.0 

-0.2 

-0.2 

B-8 

21164.7 

-4.7 

0.0 

0.0 

0.1 

0.1 

-0.0 

-0.4 

B-9 

23710.3 

-0.3 

-0.0 

-0.0 

0.1 

0.0 

-0.1 

-0.1 

B-10 

28801.4 

-1.4 

0.1 

0.1 

-0.9 

0.1 

-0.1 

-0.1 


Table 1. Timing accuracy of the background subtraction methods for the thirty simulated transients. The transients are grouped by type rather than 
chronological order, with I, E, B referring to ingreses, egresses, and bursts. The columns A 4 through Atf correspond to the error in the timings: 
b= constant cadence, c=direct subtraction, d=weighted blocks, e=weighted cells, f=weighted photons, g=Iterated Bayesian Blocks method with 
original weighting, h=Iterated Bayesian Blocks method with with alternate weighting. 


It would be desirable to have a method for dealing with 
blocks of negative count rate that has a better theoretical jus- 
tification. Sim ilar questions have previously been considered. 
iLoredol (ll992lPl gives a posterior probability distribution for the 
source count rate, independent of the count rate of the back¬ 
ground, which could potentially be used to ensure p ositive coun t 
rates in all potential blocks. A method developed bv IZechl (Il989h 
gives the probability of the source having a positive count rate s 
given the number of total (source -H background) photons de¬ 
tected, and the background count rate. This can be used to esti¬ 
mate the largest plausible source count rate. However, it is not 
immediately obvious how to incorporate either o f these proce¬ 
dures into the fitness function of Equation 19 in IScargle et alJ 
(l2013h . Both procedures also require summations over all the 
photons in a block, increasing the running time of the algorithm 
from 0(n^) to 0(rP) or worse. 

4.8. Background exposure area 

We investigated the effect of changing the exposure area of the 
background by repeating the “weighted photons” method on two 
new simulated observations, using background extraction areas 
of 2.062 and 8.246 as large as the source extraction region, that 

' A longer and more comprehensive ver¬ 
sion of this book chapter is available at 
http: //WWW. astro. Cornell. edu/staf f/loredo/bayes/t j 1. html 


is, half and double the exposure area we have been using up to 
now. The background photons were weighted according to these 
new exposure areas. 

For the three background exposure areas, we found that the 
smallest one produced 1,302 change points, the middle one pro¬ 
duced 971 change points, and the largest background exposure 
area produced 742. The accuracy in timing the thirty transient 
events was not adversely affected by increasing or decreasing 
the background exposure area, but the larger this area, the less 
spurious change points were produced. 

5. DISCUSSION AND CONCLUSION 

We have investigated the ability of the Bayesian Blocks algo¬ 
rithm to accurately and precisely time transient events. With 
an adjustment that moves the change point toward the block of 
higher count rate, the algorithm is able to determine the locations 
of instantaneous increases or decreases in count rate with essen¬ 
tially no systematic bias and with uncertainties similar in size to 
the interval between individual photons. 

When the change in count rate is not instantaneous, we have 
characterised the ability of the algorithm to resolve the period of 
varying count rate and, if it cannot be resolved, found that the 
change point is placed near the block of lower count rate. These 
observations have important implications to applications such as 
measuring the times and durations of eclipses. 
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We have incorporated background subtraction into the 
Bayesian Blocks algorithm in several different ways, and tested 
the alternatives against simulated source and background signals 
containing bursts and eclipses of known position. Even when the 
observation is dominated by large, rapidly varying, background 
contamination it is possible to recover the transient times with 
excellent accuracy. In Table[T]we show the time locations and na¬ 
ture (e.g., 1-3 is the third eclipse ingress) of the transients, as well 
as the error (the distance to the nearest bin edge found by each of 
the background subtraction methods). The results are also shown 
graphically in Figure [T3] Note that, because we know before¬ 
hand where the transients are, we can unambiguously identify 
the nearest bin edge. In real applications it may be difficult to 
determine which bin edge actually marks the beginning of the 
transient, so it is likely that methods producing many needless 
bin edges or change points will actually perform worse on real 
data than in this paper. Visual inspection of the results is always 
called for. 



Fig. 13. Timing errors of the thirty transient events, for the different 
background subtraction methods. The letter identifiers are the same as 
in Table [T] and Figure [12] b) constant cadence, c) direct subtraction, 
d) weighted blocks, e) weighted cells, f) weighted photons, g) iterated 
Bayesian Blocks with original weighting, h) iterated Bayesian Blocks 
with alternative weighting. The transient events have been divided by 
type. The constant cadence method finds all types of transients with 
precision limited by the bin size. All other methods perform well for 
bursts, but only the “weighted photons” (f) approach is reliably accurate 
at timing eclipse ingresses and egresses. The otfset points in row (f) 
show the timing errors for this method if the locations of the change 
points are not adjusted according to Equation 0 Dotted vertical lines 
indicate a timing error of zero. 

The “weighted photons” approach ( ^4.5b is clearly best for 
background subtraction. As shown in Table [1] it correctly recov¬ 
ers all transients to within a few seconds and most to within half 
of a second. However, it is prone to producing occasional blocks 
of very short lengths and implausibly high count rates (see Fig¬ 
ures [T0| and \V2\ . but if this drawback can be tolerated then the 


“weighted photons” adaption to the Bayesian Blocks algorithm 
is well suited for transient timing. All alternatives perform well 
in locating the bursts, but only the “weighted photons” approach 
reliably recovers the eclipse ingresses and egresses. The rea¬ 
son for this can be seen in Figure [T2j a pulse of background 
contamination coincides with an eclipse ingress, and this causes 
the “weighted’ blocks” and unmodified iterated Bayesian Blocks 
methods to miss the ingress entirely, and most of the other alter¬ 
natives can recover it only approximately. 

The effect of adjusting the change points according to Equa¬ 
tion |3] is also shown in Figure [T3] The offset points for the 
“weighted photons” row show the effect of not performing this 
adjustment. There is little visible difference, but the timings of 
some of the ingresses are visibly improved by performing the 
change point adjustment. 

We included a slowly rising and falling background in the 
simulated observation to determine if this causes a systematic 
offset in the timings but, from Table (Tj there is no indication of 
such an effect in any of the alternatives we considered. 

The modifications to the Bayesian Blocks algorithms devel¬ 
oped in this paper can be generalised to deal with simultane¬ 
ous observations with different instruments, each with their own 
source and background extraction regions. One simply assigns 
every photon from all photon series a weight according to the 
instrument’s extraction area and sensitivity. 

The unmodified Bayesian Blocks algorithm is invariant un¬ 
der choice of time unit. That is, it makes no difference to the 
number and placement of the change points if the photon ar¬ 
rival times are expressed in units of seconds, hours, orbital phase, 
etc. We have found that this useful property is not preserved for 
modifications to the Bayesian Blocks algorithm in which cells 
or blocks can have negative weight, probably due to the way 
we avoid taking the logarithm of a negative photon count. The 
higher the count rate, the more spurious change points will be 
produced. If the unit of time is scaled to achieve typical count 
rates of 1-10, the results should still be acceptable. 
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